Beyond power-law density scaling: Theory, Simulation, and Experiment 



O 

(N 

O 

Q 



o 



I 

O 

o 



> 

(N 
O 



X 



Lasse B0hling*, Trond S. Ingebrigtsen* , A. Grzybowski , M. Paluch", Jeppe C. Dyre*, and Thomas B. Schr0der* 
* DNRF Centre "Glass and Time", IMFUFA, Department of Sciences, 
Roskilde University, Postbox 260, DK-4000 Roskilde, Denmark, and 
"Institute of Physics, University of Silesia, ul. Uniwersytecka 4, 40-007 Katowice, Poland. 

(Dated: December 8, 2011) 

Supercooled liquids are characterized by relaxation times that increase dramatically by cooling 
or compression. Many liquids have been shown to obey power-law density scaling, according to 
which the relaxation time is a function of density to some power over temperature. We show that 
power-law density scaling breaks down for larger density variations than usually studied. This is 
demonstrated by simulations of the Kob- Andersen binary Lennard- Jones mixture and two molecular 
models, as well as by experimental results for two van der Waals liquids. A more general form of 
density scaling is derived, which is consistent with results for all the systems studied. An analytical 
expression for the scaling function for liquids of particles interacting via generalized Lennard- Jones 
potentials is derived and shown to agree very well with simulations. This effectively reduces the 
problem of understanding the viscous slowing down from being a quest for a function of two variables 
to a search for a single-variable function. 



The relaxation time of a supercooled liquid increases 
markedly upon cooling, in some cases by a factor of ten or 
more when temperature decreases by just one percent. [J- 
[Tl| This phenomenon lies behind glass formation, which 
takes place when a liquid upon cooling is no longer able to 
equilibrate on laboratory time scales due to its extremely 
long relaxation time. It has long been known that in- 
creasing pressure at constant temperature increases the 
relaxation time in much the same way as does cooling 
at ambient pressure. Only during the last decade, how- 
ever, have large amounts of data become available on how 
the relaxation time varies with temperature and density. 
Following pioneering works by T611e,[l2j it was demon- 
strated by Dreyfus et aL,[l3| Alba-Simionesco et aL.flty. 
as well as Casalini and Roland, [l5[ that for many liq- 
uids and polymers the relaxation time is a function of a 
single variable. Roland et al. 16.] reviewed the field, and 
demonstrated that for a large number of molecular liquids 
and polymers the relaxation time is a function of p 1 /T ', 
where 7 is an empirical material-dependent parameter. 
For recent works on this "power-law density scaling" or 
"thermodynamic scaling" see, e.g., Refs. ll7H20l The con- 
sensus is now that van der Waals liquids and most poly- 
mers conform to the scaling, whereas hydrogen-bonding 
liquids disobey it. 

A standard model in simulation studies of vis- 
cous liquids is the Kob-Andersen binary Lennard- Jones 
(KABLJ) mixture. plj which can be cooled to a highly 
viscous state and only crystallizes for extraordinarily 
long simulations. [22| The system consists of 80% large 
Lennard-Jones (LJ) particles (A) interacting strongly 
with 20% smaller LJ particles (B).[H[ The KABLJ mix- 
ture was shown by Coslovich and Roland[24j to obey 
power-law density scaling to a good approximation with 
7 = 5.1 for the density range p = N/V — 1.15 to 
p = 1.35, whereas Pedersen et aL[25| used the slightly dif- 
ferent value 7 = 5.16 to scale the density range 1.1 to 1.4. 
Figure [1] demonstrates, however, that power-law density 
scaling breaks down when considering a larger density 
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FIG. 1. Breakdown of power-law density scaling for 
the reduced structural relaxation time f a in the KABLJ 
mixture, f a = p 1 ^ 3 (kBT/m) 1 ^ 2 T a , where r a is the time 
at which the self-intermediate scattering function (F s (q,t), 
q = 7.25(p/1.2) x ' 3 ) for the A particles has decayed to 1/e. 
Molecular dynamics (MD) simulations in the NVT ensemble 
(N = 1000) were performed using RUMD, a MD package op- 
timized for state of the art GPU computing. [26| f a is plotted 
for three isochores as a function of the density-scaling variable 
p 7 /T, where 7 is an empirical scaling parameter. Left panel: 
7 = 4.90 collapses the data for the two lowest densities. Right 
panel: 7 = 4.45 collapses the data for the two highest densi- 
ties. It is not possible to find a single exponent that collapses 
all the data; even though the two exponents differ by only 
10%, power-law scaling with a single exponent clearly fails. 



range. Relaxation time data for the isochores p = 1.2 
and p = 1.6 collapse very well using 7 = 4.90, whereas 
the isochores p — 1.6 and p — 2.0 collapse using 7 = 4.45; 
in both cases the third isochore deviates significantly. 
In the following we show that power-law density scaling 
is an approximation to a more general form of scaling, 
which is derived from the theory of isomorphs. p7l [2q 
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We further show that given that the two lowest densities 
of Fig. 1 obey power-law density scaling with 7 = 4.90, 
the isomorph theory predicts that the two highest densi- 
ties scales with 7 = 4.45, as seen in Fig. [TJ 

What causes power-law density scaling and its break- 
down for large density variations? A justification of den- 
sity scaling may be given by reference to inverse power- 
law (IPL) potentials (oc r~ n ) where r is the distance 
between particles. For such unrealistic, purely repulsive 
systems density scaling is rigorously obeyed with 7 = n/3 
29]. Assuming that power- law density scaling reflects 
some sort of underlying effective power-law potential, the 
scaling exponent 7 can be found from the NVT equilib- 
rium fluctuations in the potential energy U and the virial 
W = pV - NkT (IPL potentials have W = (n/3)U) as 
follows: 

1 ((AC/) 2 ) 1 ) 

This was confirmed for the KABLJ mixture by Refs. [13 
and [H. Ref. further supported this "hidden scale 
invariance" explanation by demonstrating that for the 
investigated density range the dynamics and structure 
of the KABLJ mixture is accurately reproduced by an 
IPL mixture with exponent chosen in accordance with 
Eq. Q. 

The hidden scale invariance is not just a feature of 
the KABLJ mixture but of "strongly correlating liq- 
uids" in general. [30l - [33l | These are defined by having 
strong correlations between equilibrium NVT fluctua- 
tions of the potential energy and the virial (correla- 
tion coefficients larger than 0.9). Also molecular models 
can be strongly correlating; examples include the Lewis- 
Wahnstrom model of ortho-terphenyl and an asymmet- 
ric dumbell model. Both models are strongly corre- 
lating and obey power-law density scaling with expo- 
nents consistent with Eq. (JTJ for density increases of 8% 
and 16% respectively. [34j Very recently Gundermann et 
al. [20J investigated the van der Waals liquid tetramethyl- 
tetraphenyl-trisiloxane, and gave the first experimental 
confirmation of the relation between the power-law scal- 
ing exponent and Eq. ([1]). 

For any potential that is not an IPL potential the ex- 
ponent 7 as calculated from Eq. (fTJ depends on the state 
point. Power-law density scaling corresponds to disre- 
garding this state-point dependence. It is thus not sur- 
prising that power-law density scaling breaks down when 
large density changes are involved, but interestingly a 
more fundamental and robust form of scaling can be de- 
rived as we now proceed to show. 

Strongly correlating liquids have "isomorphs" in their 
phase diagram, which are curves along which structure 
and dynamics in reduced units are invariant to a good 
approximation. [13, HI] The invariance of structure im- 
plies invariance of the configurational ( "excess" ) entropy, 
S ex , thus the isomorphs are nothing but the configura- 
tional adiabats. Ref. [13 discussed in detail the conse- 
quences of a liquid having isomorphs and also showed 



that a liquid is strongly correlating if and only if it 
has isomorphs to a good approximation. The prec ise 
statistical- mechanical definition of an isomorph [27] is 
that this is an equivalence class of state points, where 
two state points are termed equivalent ( "isomorphic" ) if 
all pairs of physically relevant micro-configurations of the 
two state points, which trivially scale into one another, 
have proportional configurational Boltzmann's factors. 
From this single assumption several predictions can be 
derived, including isomorph invariance of structure and 
dynamics in reduced units and that jumps between iso- 
morphic state points take the system instantaneously to 
equilibrium. |27| 

Letting R denote a micro-configuration (all particle co- 
ordinates) of a reference state point (p*,T*), the condi- 
tion for state point (p, T) to be isomorphic with (p*, T*), 
i.e., the proportionality of pairs of Boltzmann's factors, 
can, by taking the logarithm and rearranging, be ex- 
pressed as: 

u(jrin) = ^rU(R) + K, p = P / P * (2) 

where if is a constant that only depends on the two 
state points. Equation ([2]) is the basis of the so called 
direct isomorph check: [271 ] a) draw micro-configurations 
R from a simulation at (p*,T*), b) evaluate the potential 
energies of these configurations scaled to density p, and 
plot them in a scatter-plot against the potential energies 
at p* . If a state point (p, T) exists that to a good ap- 
proximation is isomorphic with (p*, T*), this scatter plot 
will be close to a straight line and T can be found from 
the slope. 

In the following we consider systems where the inter- 
action potential between particles i and j is given by a 
sum of inverse power laws: 

=^2 e n,ij ( — ) ■ (3) 

This includes the standard 12-6 LJ potential, but also, 
e.g., potentials with more than two power-law terms. We 
note that some systems interacting via Eq. ([3]) will not be 
strongly correlating and thus not have good isomorphs. 
In the following properties are derived for those systems 
that do have good isomorphs. 

The total potential energy of a given micro- 
configuration R at density p* is a sum over contribu- 
tions from the power-law terms, U — ^2 n U n - When 
scaling R to the density p, keeping the structure invari- 
ant in reduced units, each power-law term is scaled by 
pa = (p/p*)"/ 3 , and the potential energy at the new 
density U' = U (p- x / 3 R) is:(H 

U' = J2p f U n . (4) 

n 

Thus the linear regression slope of the U', [/-scatter plot 
in the direct isomorph check is given by (where all aver- 



3 



ages refer to the reference state point (p*,T*)): 
(AU'AU) 
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Using Einstein's fluctuation formula for the excess iso- 
choric heat capacity and the corresponding formula for 
the "partial" heat capacities (which can be negative), 
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we get an expression for the new temperature T relative 
to the reference temperature T* (compare Eq. @): 
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Since = J2 n ^Tn the number of parameters in the 
scaling function g(p) is one less than the number of terms 
in the potential (Eq. ©). In particular, for the standard 
12-6 LJ potential, g(p) contains just a single parameter: 
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Using that U 12 = W/2 - U for 12-6 LJ systems, [H] g{p) 
can conveniently be expressed in terms of 7* defined as 
Eq. ([!]) evaluated at the reference density p*: 



g{p) 



p 4 ( 7 */2 -l)-p 2 ( 7 */2 - 2) . 



(9) 



g(p) provides a new and convenient method for numer- 
ically tracing out an isomorph - previously this could 
only be done by changing density by a small amount, 
e.g., 1%, and then adjusting temperature to keep the ex- 
cess entropy constant, using that 7 (Eq. (JT}) also can be 
expressed as:[27| 
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d\nT 
dlrip 



(10) 



It is a prediction of the isomorph theory that 7 depends 
only on density. [27J, |35[ This means that the same differ- 
ential equation, Eq. (fl0|) . determines the temperature on 
all isomorphs, implying that g(p) is the same for all iso- 
morphs - what changes between different isomorphs is . 
Thus g(p)/T is an isomorph invariant (compare Eq. ([7])). 
which can be used as a scaling variable for the reduced 
relaxation time f that is also an isomorph invariant f27| 
f = f(g(p)/T). This form of scaling was first proposed 
by Alba-Simionesco et aZ.[lJ| Here a theoretical deriva- 
tion has been provided, as well as an explicit expression 
for g(p) for systems interacting via generalized LJ poten- 
tials (Eq. (|5J|). We note further that since the solid-liquid 
coexistence lines of st rongly correlating liquids are pre- 
dicted to be isomorphs [2 71 128|| . this immediately explains 
the recent observation of Khrapak and coworkers that 
the solid-liquid coexistence line of an LJ liquid is given 
by (Cip 4 - C 2 p 2 ) /T = Const. [H 
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FIG. 2. Four different isomorphs in the KABLJ mixture, 
each generated from the condition g(p)/T — Const. Densities 
range from 1.2 to 2.0, and g(p) = p 4 (7*/2 - 1) -p 2 (7»/2 - 2) 
(Eq. (|9}) with 7* = 4.59 determined from the scaling in Fig.[T] 
see text (p = p/p*, p* = 1.6). (a) Self part of intermediate 
scattering functions in reduced units, (b) Mean-square dis- 
placements in reduced units. The data collapse confirms that 
true isomorphs have been identified. 



From the scaling f = f(g(p)/T) follows that pairs of 
isochores obey power-law density scaling, since a power- 
law Ap 1 that is equal to g(p) at the two densities in- 
volved always exists. This is indeed what is seen in 
Fig. [U Choosing p* = 1.6 as reference density, the scal- 
ing in the left panel of Fig. Q] corresponds to 5(1.2/1.6) = 
(1.2/1.6) 4 90 which via Eq. © leads to 7* = 4.59. Us- 
ing this value we find #(2.0/1.6) = 2.70 = (2.0/1. 6) 4 45 . 
Thus from one scaling exponent in Fig. [T] the other is 
uniquely predicted. Moreover, the value 7, = 4.59 is 
consistent with what is found by evaluating Eq. (fTJ) at 
the reference isochore p* = 1.6 in the temperature range 
T = 1.7 to 5, which leads to values of 7* decreasing from 
4.6 to 4.5. In the following figures reporting results for 
the KABLJ mixture, we use 7* = 4.59 as estimated from 
the left panel of Fig. [TJ i.e., no further fitting or adjust- 
ment of parameters was applied. 

As mentioned, the scaling function g(p) was derived 
assuming that good isomorphs exist. In Fig. [5] we test 
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FIG. 3. Reduced relaxation times for the KABLJ mixture 
scaled according to the isomorph theory (same data as in 
Fig. [TJ. The scaling function g(p) is the same as in Fig. [5] 
(Eq. 7» = 4.59). Inset: Comparing 7 computed from 
simulations (Eq. fTJ) to the prediction of the isomorph the- 
ory, 7 = dlng/dlnp (black curve). Isomorphs denoted by 
a reduced relaxation time are each generated from a single 
reference point using Eqs. (0 and (J5J, with '±' quantifying 
the resulting variation of the reduced relaxation time on the 
isomorph. Isochores are plotted with the same symbols as in 
the main figure. 



this for the KABLJ mixture using the most sensitive iso- 
morph invariant - the dynamics of viscous states. State 
points with the same g(p) /T, predicted to be on the same 
isomorph, are seen to have very similar dynamics even 
though density varies from 1.2 to 2.0. 

Figure[3]tests the proposed scaling for the KABLJ mix- 
ture using the data of Fig. [TJ Clearly the new form of 
scaling works well. Combining Eq. (|10p with the defi- 
nition of g(p) (Eq. (J7J) shows that 7 is the logarithmic 
derivative of g(p), 7 = dlng/dlnp. The inset of Fig. [3] 
demonstrates that this prediction agrees well with sim- 
ulations, even when going to large densities where the 



purely repulsive 



limit is approached. 



We now turn briefly to model molecular systems. In 
this case it is still a prediction of the isomorph theory 
that an expression of the form g(p)/T is the right scal- 
ing variable [35j], but we do not have an explicit expres- 
sion for g(p). Figure U demonstrates how power-law den- 
sity scaling breaks down for the Lewis- Wahnstrom model 
of ortho-terphenyl and an asymmetric dumbell model 
when considering larger density changes than previously 
studied. [HI Like in Fig.[TJ power- law scaling works when 
considering pairs of isochores, consistent with the right 
scaling variable being of the form g(p)/T. The insets of 
Fig. [4] test the isomorph prediction that 7 to a good ap- 
proximation is a function of density only, the assumption 
used to derive the new scaling. The prediction agrees well 
with simulations: 7 is found to be much more dependent 
on density than on temperature. For more results on 
isomorphs in these model molecular liquids see Ref. H3- 

Finally, we present experimental results for the two van 
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FIG. 4. Breakdown of power-law density scaling in two 
molecular models. In accordance with the scaling derived 
in the present work, power-law scaling does work when 
considering only pairs of isochores. a) Asymmetric dumb- 
ell model, b) Lewis- Wahnstrom model of ortho-terphenyl 
(OTP). Insets: (dln-y/dlnT) p plotted against InT (circles) 
and (din y/d In p)r plotted against lnp (stars). Consistent 
with the isomorph theory, 7 is found to be much more depen- 
dent on density than on temperature. 



der Waals liquids dibutyl phthalate (DBP) and decahy- 
droisoquinoline (DHIQ), using larger density increases 
than usually studied in scaling experiments (20% and 
18% respectively). Using earlier reported dielectric data 
for DBP[H] and DHIQ [39] and values of the Tait equa- 
tion of state parameters for DBP[i(| and DHIQ[4l|, we 
find that the isochronal dependences log 10 T vs log 10 p 
determined at given structural relaxation times in re- 
duced units, t — TVm^(kBT/m) 1 / 2 1 where v m and m 
are molecular volume and mass, are nonlinear (Fig. [SJb) 
and (e)). This implies breakdown of power-law density 
scaling. The isochrones can be superposed after vertical 
shifting, however (Fig. [Sic) and (f)), which implies that a 
scaling variable exists of the form g(p) /T . The isochrones 
can be described by a phenomenological form of the seal- 
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FIG. 5. Deviation from power-law scaling in DBP (a-c) and DHIQ (d-f). (a) and (d) Density dependence of isobaric and 
isothermal structural relaxation times r~in reduced units. Solid lines represents separate fits to the modified Avramov model. [43] 

(b) and (e) The isochronal dependences log 10 T vs log 10 p determined at given f. Fits are done to all isochrones simultaneously. 

(c) and (f) The isochrones vertically shifted by the fitted values A(f). The deviations from power- law scaling (straight dashed 
lines) are evident. The fitted g(p)'s corresponds to 7 values increasing from 2.6 to 3.9 (DBP) and from 2.0 to 4.3 (DHIQ). 



ing function log 10 g(p) = A\ log 10 p + A 2 log 10 p, chosen 
simply to take into account the curvature. In the case 
of DBP, our results are consistent with those reported 
by Niss et al.JH] We conclude that for DBP and DHIQ 
power-law density scaling breaks down at large density 
variations in the way predicted by the isomorph the- 
ory. Interestingly, the density dependence of 7 is stronger 
than for the model systems, and in the opposite direction; 
for the experimental systems 7 increases with density. 

What are the perspectives of our findings? Based 
on the theory of isomorphs in dense liquids we have 
now a form of density scaling that is more fundamen- 
tal and more robust than power-law density scaling, and 
which is consistent with both simulations and experi- 
ments. This "isomorph scaling" - originally proposed 
by Alba-Simionesco et al. 1J] - is predicted to apply for 



all strongly correlating liquids, i.e., van der Waals and 
metallic liquids, but not, e.g., for hydrogen-bonding liq- 
uids. Our results should not be used to abandon power- 
law density scaling - it is a useful approximation to iso- 
morph scaling when the scaling function g(p) is unkown 
and only moderate density changes are considered. Un- 
der these conditions the isomorph theory predicts that 
power-law density scaling works with an exponent deter- 
mined by Eq. ([1}. Isomorph scaling provides a deeper 
understanding of why - and when - power-law density 
scaling works. 

Isomorph scaling has important consequenses regard- 
ing the most fundamental open question in the field of 
viscous liquids and the glass transition: what controls the 
dramatic viscous slowing down? In general this question 
has to be considered as a search for a physically justified 
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function of two variables, temperature and density (or 
temperature and pressure). Our results imply that this 
problem is now effectively reduced to a search for a func- 
tion of a single variable, at least for the class of strongly 
correlating liquids. This is particularly striking for LJ 
type systems like the KABLJ mixture, where we have a 
prediction for the scaling function, that agrees very well 
with simulations for much larger density variations than 
usually considered. The fact that the LJ scaling func- 
tion contains just a single parameter - i.e., no more pa- 
rameters than power-law density scaling - confirms that 
isomorph scaling is more fundamental and not merely a 
higher-order approximation compared to power-law den- 



sity scaling. Isomorph scaling must be taken into account 
in theories of the viscous slowing down: since the relax- 
ation time in reduced units obeys isomorph scaling, any 
quantity proposed to control the relaxation time must 
also obey isomorph scaling. 
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